#include "PSD.h"
#include "mainwindow.h"
#include<QDebug>

PSD::PSD(int n,int fs,QString filename)
{
    N = n;
    Fs = fs;
    Filename = filename;
}

PSD::~PSD()
{
    return;
}

QVector<QVector<double>> PSD::seepsd()
{
    QFile file(Filename);
    QVector<double> array; //用来保存信号
    QVector<double> array_x; //画图的时候使用，array_x作为横坐标，array作为纵坐标
    if(file.open(QIODevice::ReadOnly))
    {
        QTextStream stream(&file);
        while(!file.atEnd())
        {
            double buf;
            QStringList list=stream.readAll().split("\n");
            QListIterator<QString> li(list);
            while(li.hasNext())
            {
                buf=li.next().toDouble();
                double temp = buf ;
                array.append(temp);
            }
        }
    }
    //进行自相关变化
    int len = array.size();
    QVector<double> array_corr;//将信号进行自相关转变
   /* for(int i=0;i<=2*len+1;i++)
    {
        array_corr.append(0);
    }
    for(int i=0;i<=len;i++)
    {
        for(int j=0;j<=i;j++)
        {
            array_corr[i]+=array[j]*array[len-1-i+j];
        }
    }
    for (int i=len;i<(len<<1)-1;i++)
    {
        array_corr[i]=array_corr[(len<<1)-2-i];
    }
    */
    //进行傅里叶变换
    fftw_complex *in, *out;
    in = (fftw_complex *)fftw_malloc(sizeof(fftw_complex) * N);
    out = (fftw_complex *)fftw_malloc(sizeof(fftw_complex) * N);

    for (int i = 0; i < N; i++)
    {
        in[i][0] = array[i];
        in[i][1] = 0;
        array_x.append(i);
    }

    fftw_plan my_plan;
    my_plan = fftw_plan_dft_1d(N, in, out, FFTW_FORWARD, FFTW_ESTIMATE);
    fftw_execute(my_plan);

    QVector<double> data_power;
    QVector<double> data_power_x;
    QVector<double> data_freq;


    //功率谱密度估计
    for( int i = 50; i < 100000; i++) //这里的i的区间可以自由选择。
    {
        double tempx = (out[i][0] * out[i][0] + out[i][1] * out[i][1])/(N);
        tempx = sqrt(tempx);
        data_power.append(tempx);
//      data_cali_psd_frequency.append((2*Pi*i)/(N*T));
        data_freq.append(i*6.28*48/100000);
        data_power_x.append(i-50);
    }
    fftw_destroy_plan(my_plan);
    fftw_free(in);
    fftw_free(out);
    QVector<QVector<double>> data;
    data.append(data_freq);
    data.append(data_power);
    return data;
}
